Wavepacket dynamics in energy space of a chaotic trimeric Bose-Hubbard system 
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I. INTRODUCTION 

Understanding the intricate behavior of bosonic many- 
body systems has been a major challenge for leading re- 
search groups over the last years. Without doubt, the 
theoretical interest was strongly enhanced by recent ex- 
perimental achievements in handling (ultra-)cold quan- 
tum gases: Namely, since the celebrated realization of 
atomic Bose- Einstein condensates (BEG) in periodic op- 
tical lattices (OL) [TJ El |3j H] and the creation of "atom 
chips" |S] [SJ 13 [H] we have versatile tools at hand which al- 
low for an unprecedented degree of precision as far as ma- 
nipulation and measurement of the atomic cloud is con- 
cerned. While this has led on the one hand to novel, con- 
crete applications of quantum mechanics like e.g. atom 
interferometers P [TOJ HI] and lasers [2 H [H [H [H] it 
also enabled us to investigate complex solid-state phe- 
nomena, such as the Mott-Insulator to superfluid transi- 
tion [TS] or the Josephson effect [3]. 

Beside these advances, our understanding of bosonic 
many-body systems is still very limited once we consider 
an (external) driving: In the framework of BECs in OLs 
this can be, e.g., a modulation of the potential height or 
a tilting of the lattice. Due to the time dependence of 
the driving parameter, the energy of the system is not 
a constant of motion. On the contrary, the system ex- 
periences transitions between energy levels and therefore 
absorbs energy. This irreversible loss of energy is know 
as dissipation [ini [T71 [THl [T5] . The classical dissipation 
mechanism is by now well understood |16] while quantum 
dissipation still poses some challenges. In order to get a 
better insight into the problem the main task is to build a 
theory for the time-evolving energy distribution. Apart 
from being of fundamental interest, such a theory will 
also shed a new light on recent experiments with BECs 
amplitude-driven OLs [201 EI] ■ It has been pointed out 
that the energy absorption rate measurements can be 
used to probe the many-body excitations of the system 

EH E3 El ESI Ell- 
in the present work, we approach the problem of quan- 



tum dissipation by studying the quantum dynamics of 
interacting bosons on a ring-shaped lattice consisting of 
three sites (trimer). Specifically, we will analyze the sys- 
tem's response to a rectangular pulse of finite duration 
t that perturbs the coupling fco fc = fco + between 
adjacent sites. In the framework of OLs this corresponds 
to a sudden change in the intensity of the laser field and 
is readily achieved in the experiment [SUl Ell • The asso- 
ciated dynamical scenario known as wavepacket dynam- 
ics [m EH EH ES] is one of the most basic non-trivial 
evolution schemes. Its analysis will pave the way to un- 
derstand more demanding evolution scenarios and ulti- 
mately the response of interacting bosons under persis- 
tent driving. 

The minimal quantum model that describes interact- 
ing bosons on a lattice M wells is the Bose-Hubbard 
Hamiltonian (BHH) , which incorporates the competition 
between kinetic and interaction energy of the bosonic sys- 
tem. The BHH is based on a M-mode approximation and 
hence its validity is subject specific conditions discussed 
in Rcfs. [llSnilSIlISS (see also Section |ll|. As far as the 
BHH is concerned, the two-site system (dimcr) has been 
analyzed thoroughly from both the classical (mean-field) 
[33, iAj j55^ and the purely quantum viewpoint [33 , 36 , 371 
and many exciting results were found including their ex- 
perimental realization |38j . 

As a matter of fact, the dimer is integrable since the 
BHH has two conserved quantities, the energy and the 
number of bosons. The addition of a further site -yielding 
cither a linear chain (open be) or a ring (periodic bc)- 
is sufficient to make the resulting system (trimer) non- 
integrable and thus leads to (classically) chaotic behav- 
ior. Here we consider a three-site ring |87j which can 
be experimentally realized using optical lattice or micro 
trap technology [6l III HI |39l . For example, an optical po- 
tential in a ring configuration can be achieved by letting 
a plane wave interfere with the so-called Laguerre-Gauss 
laser modes as described in [40i |41j . Another possibility 
to experimentally create a three-site ring BEG trap jMl 
is given by a combination of the methods described in 
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Refs. Il^l US]: In the experiment of (12] a trapping po- 
tential is partitioned into three sections by a central re- 
pulsive barrier created with blue-detuned laser light that 
is shaped to segment the harmonic oscillator potential 
well into three local potential minima. For better opti- 
cal resolution (up to 1.2/im), and control of the coupling 
between the three condensates one can substitute the de- 
tuned laser source of Ref. [Hj with the one used in |33j • 

The motivation to study the quantum trimcr is 
twofold: That is, while remaining simple enough to allow 
for a thorough analytical study, it displays a whole new 
class of complex behaviors which are typical for longer 
lattices consisting of many sites. The trimer has been 
studied quite extensively in the classical (mean-field) 
regime HSJ |35] . Less attention was paid to the analy- 
sis of the quantum trimer [HllTlllHlllilMJEIllSa. As 
a matter of fact, the majority of these studies is focused 
on the statistical properties of levels [47l HHl [Sj while 
recently an analysis of the shape of eigenstates was per- 
formed in Ref. [29 . However, the knowledge of spectral 
and wave function statistics is not enough if one wants to 
predict the dynamical behavior of a system. 

In our study we combine three theoretical approaches: 
On the one hand, we will use linear response theory 
(LRT) which constitutes the leading framework for the 
analysis of driven systems [TB]. On the other hand, 
we employ an improved random matrix theory (IRMT) 
modeling. Although random matrix theory (RMT) was 
proven to be a powerful tool in describing stationary 
properties (like level statistics |1H1 [S3] and eigenfunctions 
[2S]), its applicability to the description of wavepacket 
dynamics is not obvious [HI . The latter involves not 
only the knowledge of the statistical properties of the 
two quantities mentioned above but also the specific cor- 
relations between them. Finally, we will investigate the 
validity of semiclassical methods to describe the quantum 
evolution. Our analysis indicates that some moments of 
the evolving energy distribution show a remarkable level 
of quantum-classical correspondence (QCC) [2S| 
while others are strongly dominated by quantum inter- 
ference phenomena. 

The structure of this paper is as follows: in the next 
section, we introduce the Bose-Hubbard Hamiltonian 
that mathematically describes a quantum three-site ring- 
lattice. We identify its classical limit, leading to the 
discrete nonlinear Schrodinger equation and derive the 
classical equations of motion. In Section III we discuss 



the notion of wavepacket dynamics and introduce the ob- 
servables studied in the rest of the paper. We begin our 
analysis with the statistical properties of the spectrum 



and of the matrix elements of the BHH (Section IV I 



This study allows us to introduce an IRMT modeling 
which is presented in Subsection |IVD[ In Section |V] we 
extend our previous analysis on the parametric evolution 
of the eigenstates of the BHH [29] by comparing the ac- 
tual quantum mechanical calculations with the results of 
the IRMT modeling. We introduce the concept of para- 
metric regimes [29j and show how it can be applied to 



analyze the parametric evolution of the local density of 
states (LDoS) [HI [29l EU • We then turn to the dynamics 
of the BHH (Section [Vl| and extend the notion of regimes 
to the wavepacket dynamics scenario. The predictions 
of LRT, IRMT modeling and semiclassics are compared 
with the exact quantum mechanical calculations for the 
trimeric BHH model. We find that the energy spreading 
SE(t) shows a remarkable quantum-classical correspon- 
dence which is independent of the perturbation strength 
dk. In contrast, other observables are sensitive to quan- 
tum interference phenomena and reveal QCC only in the 
semiclassical regime. The latter can be identified with 
the non-perturbative limit associated with perturbations 
Sk > Skprt- Section VII summarizes our findings. 



II. THE BOSE-HUBBARD HAMILTONIAN 

The mathematical model that describes interacting 
bosons in a (three-site) lattice is the Bose-Hubbard 
Hamiltonian, which in second quantization reads 



U ^ 

H ^ -J2Mn^-^)-kY,b^b 



J' 



n ^ 1 



(1) 



Here we consider a three-site ring configuration which 
is experimentally feasible with current optical methods 
where, for example, the trapping potential is created by 
letting a plane wave interfere with the so-called Laguerre- 
Gauss laser modes as described in [ID]. The operators 
hi = SJSi count the number of bosons at site i. The anni- 
hilation and creation operators bi and obey the canon- 
ical commutation relations = Sij. In the BEC 
framework, k = ko + Sk, is the coupling strength between 
adjacent sites and can be controlled experimentally 
(in the context of optical lattices this can be achieved by 
adjusting the intensity of the laser beams that create the 
trimeric lattice), while U — ^Trfi^as/mVcS describes the 
interaction between two atoms on a single site (to is the 
atomic mass, as is the s-wave scattering length of atoms 
which can be either positive or negative, and Vcs is the 
effective mode volume). It is interesting to note that the 
BHH also appears in the context of molecular physics 
where |36} I55j k represents the electromagnetic and me- 
chanical coupling between bonds of adjacent molecules 
while U represents the anharmonic softening of the 
bonds under extension. 

The Bose-Hubbard model for M sites is based on a 
A/- mode approximation [30 (in the limit of long lattices 
this corresponds to a single (lowest) band approximation 
of the OL f?). This assumption holds provided that the 
chemical potential, the kinetic energy and the interaction 
energy are too low to excite states in the higher single- 
well modes (higher Bloch bands accordingly). Therefore, 
the lattice must be very deep [301 EU ES] inducing large 
band gaps. Furthermore, the interaction energy has to 
be smaller than the single particle ground state energy. 
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so as to not considerably modify the single particle wave- 
function. A Gaussian approximation of the wavefunction 
together with a standard harmonic trap of size lO/^m and 
a scattering length Os = 5nm indicates that the BHH 
model is valid for up to several hundred bosons per trap 
[30]. 

Hamiltonian ([T]) has two constants of motion, namely 
the energy E and the number of particles N = 
Having N = const, implies a finite Hilbert-space of di- 
mension J\f ^ {N + 2){N + l)/2 [36i 113 which can be 
further reduced by taking into account the threefold per- 
mutation symmetry of the model |29j . 

For large particle numbers N ^ 1 one can adopt a 
semiclassical approach for Hamiltonian ([ij. Formally, 
this can be seen if we define rescaled creation and anni- 
hilation operators — bt/VN. The corresponding com- 
mutators [ci, Cj] — Sij/N vanish for 3> 1 and therefore 
one can treat the rescaled operators as c-numbers. Using 
the Heisenberg relations Ci \fTi exp**^' {^pi is an an- 
gle and li is the associated action [52] ), we obtain the 
classical Hamiltonian H. 

^ - ^ - ^ E - A E y^exp^(— ) , (2) 

where C/ = NV is the rescaled on-site interaction. 

The dynamics is obtained from (|2]) using the canoni- 
cal equations dli/dt = —dTi/dipi and dipi/dt = d'H/dli. 
Here t = U-t is the rescaled time. The classical dynamics 
depends both on the scaled energy E = E/UN and the 
dimensionless parameter A = k/fj [551 liSl liTl liSl 157]. 
For A the interaction term dominates and the sys- 
tem behaves as a set of uncoupled sites (also known as 
the local-mode picture [36 ) while in the opposite limit of 
A — > oo, the kinetic term is the dominant one {normal- 
mode picture [34l |57l [58] )• In both limits the motion is 
integrable while for intermediate values of A the trimeric 
BHH ([T]) has a chaotic component [15]. We point out 
that the classical limit is approached by keeping A and 
U constant while N ^ oo [25]. This is crucial in order 
to keep the underlying classical motion unaffected. 



III. PRELIMINARY CONSIDERATIONS AND 
OBJECT OF THE STUDY 

In this paper we study the trimeric BHH model ^ as a 
control parameter, the coupling strength between lattice 
sites is changed i.e. fco ^ A;o -I- 5k. In our analysis, we 
therefore consider 

H ^Hq- 5k{t)B , (3) 
where the perturbation operator B is 

^ - E ^1^^ ' (4) 
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FIG. 1: (Color online) Scheme of the wavepacket dynam- 
ics scenario: the perturbation is a rectangular pulse of dura- 
tion ti at which the measurement is done. The function f{t) 
represents the rescaled time-dependence of the perturbation 
Sk{t) = 5k X f{t) (black hue) while the red line indicates its 
time derivative f{t). 

and the unperturbed Hamiltonian Hq is given by Eq. |l]) 
with k = ko. Quantum mechanically, we work in the 
Hq eigenbasis. In this basis Hq becomes diagonal, i.e., 
Eq — Em^Smn whcrc {Em^} are the ordered eigenvalues 
and we can write 

H = Eo - . (5) 

Throughout this work we always assume that the per- 
turbed Hamiltonian T-l(k) as well as the unperturbed 
Hamiltonian 7i(fco) generate classical dynamics of the 
same nature, i.e., that the perturbation Sk ^ k — fco is 
classically small, 5k < 5k^i (see beginning of the next 
section for the definition of 5k^i). This assures the appli- 
cability of classical linear response theory (LRT). Note, 
however, that this assumption is not sufficient to guar- 
antee the validity of quantum mechanical linear response 
theory. Our aim is to identify novel quantum mechani- 
cal effects that influence the classical LRT results as the 
perturbation 5k increases. At the same time, we address 
the implications of classically chaotic dynamics for the 
trimeric BHH, and the route to quantum-classical corre- 
spondence in the framework of wavepacket dynamics. 

For later purposes it is convenient to write the per- 
turbation as 5k{t) — Sk X f{t) where 5k controls the 
"strength of the perturbation" while f{t) is the scaled 
time dependence (note that if we had f{t) (x t, i.e. persis- 
tent driving, then 5k would be the "rate" of the driving) . 
Although our focus will be on the wavepacket dynamics 
scenario where the perturbation is a rectangular pulse of 
strength 5k and duration t -see Fig. [ij for a sketch of the 
resulting step function f{t) with k{t) = fc(0)- we expect 
that the results presented here will shed some light to 
the response of BHHs in the presence of more demand- 
ing driving scenarios. 

A. Measures of the evolving distribution Pt{n\no) 

In this subsection we discuss a number of observables 
that will allow us to quantify the response of the system 
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FIG. 2: Poincare sections of the phase space belonging to the classical trimer for N = 1 and different parameter values a) 
A = 0.005, b)A = 0.05, c) A = 2. On the y-axis we plot the action I3 while on the a;-axis the difference ifi2 — f3 (in units of tt) 
is plotted. The Poincare section corresponds to the plane ipi — ips and ipi > (p2 of the energy surface E — 0.2. 



and the spreading of the energy distribution. 

We consider an initial micro-canonical preparation de- 
scribed by an eigenstate |no) of the unperturbed Haniil- 
tonian H{k{0)). Given the driving scenario k{t), it is 
most natural to analyze the evolution of the probability 
distribution 



Ptin\no) = \{n\Uit)\no)\' 



where 



U{t) 



Texp[-- / dt'H{k{t')] 



(6) 



(7) 



is the time-ordered evolution operator and 
H[k{t)]\n[k{t)]) = En[k{t)]\n[k{t)]) . By convention 
we order the states by their energy. Hence we can regard 
Pt{n\no) as a function of r = n — no, and average over 
the initial preparation (around some classically small 
energy window), so as to get a smooth distribution 
Pt{r). 

To capture various aspects of the evolving probabil- 
ity distribution Pt{n\no) we introduce here the survival 
probability defined as 



P{t)^\{no\U{t)\no)\^^Pt{no\no), 



and the energy spreading 



SEit) - /^Pt(n|no)(i;„ 



(8) 



(9) 



which probes the tails of the evolving distribution. Yet, 
the evolution of Pt{n\no) is not completely captured by 
any of these measures: As we will see in Section [VT] the 
wavefunctions can develop a "core" which is a result of 
a non-perturbative mixing of levels [29]. We therefore 
define an operative measure that reflects the creation of 
the "core", as the width SE^^^^ which contains 50% of the 
probability: 



(5£;„„(t) = [7175% - n25%]A 



(10) 



Here, A is the mean level spacing and Ug is determined 
through the equation J2n ^t('^l'^o) = 1- 



IV. STATISTICAL PROPERTIES OF THE 
TRIMERIC BHH: SPECTRA AND BAND 
PROFILE 



The dynamical properties of the classical trimer were 
thoroughly investigated in a number of papers |441 145} 
H5] . It was found that for intermediate values of the con- 
trol parameter A, the system exhibits (predominantly) 
chaotic dynamics. Some representative Poincare sections 
(corresponding to the plane tpi = ip^ and ipi > ip2 of 
the energy surface E = 0.2 of Hamiltonian ^) of the 
phase space are reported in Fig. |2] As A decreases, 
one can clearly see the transition from integrability to 
chaotic dynamics and back to integrability. We deter- 
mine the regime of predominantly chaotic motion based 
on the nature of the phase space and the power spec- 
trum C{lu) of the classical perturbation operator (the 
latter is discussed in detail in Subsection IV B I . While 



regular motion results in isolated peaks in C{lu), a con- 
tinuous (but possibly structured) power spectrum indi- 
cates chaoticity. Accordingly, the classical smallness con- 
dition Sk <^ Skci can be operatively deflned as the per- 
turbation strength that leaves C{uj) unaffected. We have 
found that for 0.04 < A = k/U < 0.2 and an energy 
interval H w 0.26 ± 0.02 the motion is predominantly 
chaotic. Choosing our parameter values to be feg = 15 
and U = 280 we find Sk^i « 20. 

In the following we will concentrate on the above men- 
tioned range of A values for which chaotic dynamics is 
observed. The main question we will address is: What 
are the signatures of classical chaos in various statisti- 
cal quantities upon quantization? As we shall see in 
the following subsections, chaos manifests itself mainly in 
two quantities; the spectral statistics of the eigenvalues 
{Em^} and the averaged profile (|Bm„p) of the pertur- 
bation operator. While the statistical properties of the 
levels have attracted some attention in the past |47l [59] , 
the traces of chaotic dynamics in the shape of the pertur- 
bation operator (|Bto„P) and the statistical properties of 
its matrix elements were left unexplored. In the next sub- 
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h=k/U 

FIG. 3: (Color online) Parametric evolution of the eigen- 
values En'' as a function of the parameter A. The number 
of bosons is A'' = 40 and the effective interaction strength is 
U = 280. In the main figure the entire spectrum is plotted 
while the inset is a magnification of the small box. One ob- 
serves a qualitative change in the spectrum as A is changed. 
See text for details. 



sections we will address these issues in detail and propose 
an improved random matrix theory modeling which takes 
our statistical findings into consideration. 



A. Energy levels 

In Fig.lSlwe plot the parametric evolution of the eigen- 
values En as a function of A for fixed effective interac- 
tion strength U = 280. From Fig. [sjone observes that 
the spectrum becomes rather regular for very large A. In- 
deed, for A — *■ oo a transformation to the normal modes 
of the system diagonalizes the Hamiltonian and yields 
an equidistant spacing of the eigenvalues [47]. In the 
local- mode limit, i.e. A — > 0, the eigenvalues of Hq are 
obtained immediately from ([T]) and are partly degenerate 
[90| . However, in an intermediate A-regime one observes 
a different behavior, namely irregular evolution and level 
repulsion (see inset) . This is a manifestation of the clas- 
sically chaotic behavior ^60| ^61^ . 

In order to establish this statement we turn to the 
statistical properties of the spectra. In particular, we will 
study the level spacing distribution V{S) [iTi i48i i53i. i59 
where 



En+1 — En 



(11) 



are the spacings of two consecutive energy levels which 
are unfolded with respect to the local mean level spacing 
A. The level spacing distribution represents one of the 
most popular measures used in quantum chaos studies 
[60|, [6T| . It turns out that the sub-?i statistical features 
of the energy spectrum of chaotic systems are "univer- 
sal" , and obey the RMT predictions [521 [53] • In contrast. 



non-universal, i.e. system specific, features are reflected 
only in the large scale properties of the spectrum and 
constitute the fingerprints of the underlying classical dy- 
namics. 

The mean level spacing A can be estimated from the 
fact that J\f cx iV^ levels span an energy window AE (x 
UN X E, around some specific energy E (see Eq. ([l])). 
Our considerations indicate the scaling relation 



A « 1.5 



U_ 

N ' 



(12) 



where the proportionality factor was found by a direct fit 
of our spectral data in the energy window around E — 
0.26 [29;. 

For chaotic systems the level spacing distribution 'P{S) 
follows the so-called Wigner surmise [BTJ IM] 



Se- 



as) 



indicating that there is a linear repulsion between nearby 
levels. Instead, for generic integrable systems there is no 
correlation between the eigenvalues and the distribution 
'P{S) is Poissonian 



rus) 



(14) 



In Fig. [4] we report some representative V{S) for lev- 
els in the energy window around E = 0.26 [9T]. One 
observes a qualitative change in the shape of 7^(5*) from 
Poissonian-like associated with very small and large A 
values to Wigner-like for intermediate values of A. 

In order to quantify the degree of level repulsion (and 
thus of chaoticity), various phenomenological formulas 
for 'P{S) have been suggested that interpolate between 
the two limiting cases (13] 14 1 (see for example [551155] ). 
Here we use the so-called Brody distribution [55] given 
by the following expression 



Vq{S)^aS'' e~^^' 



(15) 



where a = {1 + q)f3, (3 = ri+«[(2 + q)/{l + q)] and F 
is the Gamma function. The two parameters a, /3 are 
determined by the condition that the distribution is nor- 
malized with a mean equal to one [57*. The so-called 
Brody parameter q is then obtained from direct fitting of 
VqiS) to the numerically evaluated level spacing distribu- 
tion. One readily verifies that for q = 0, the distribution 
'Pq{S) is Poissonian ( 14 1 while for q — 1 it takes the form 
of ([Tl. 



The fitted values of the Brody parameter q for various 
A's are summarized in Fig. |5] We see that for very small 
and very large A the Brody parameter is small indicating 
classically regular motion while for intermediate values 
0.04 < A < 0.2 we find q ^ 1 corresponding to classically 
chaotic motion. This result is in perfect agreement with 
the predictions of the classical analysis. 
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FIG. 4: (Color online) The level spacing distribution 1^(3) of the BHH trimer for three representative values of the dimensionless 
ratio A = k/U which controls the underlying classical dynamics: a) A = 0.025 (fc = 7), b) A = 0.05 (fc = 14.5), and c) 
A = 0.35 (fc = 100). The red dash-dotted line corresponds to the Pois sonian distribution (14 1 which is expected for integrable 



systems, the solid blue line corresponds to the Wigner surmise ( 13 1 (chaotic systems) while the solid green line represents 
the fitted Brody distribution (15 1. In Figure [5] we report the fitted Brody parameter q for various values of A. The System 
corresponds to A'^ = 230 bosons and U — 280. The histograms include the 400 relevant levels around E — 0.26. 




FIG. 5: The Brody parameter q for the BHH plotted against 
the dimensionless ratio A which controls the underlying classi- 
cal dynamics. The values of q are obtained from fits to Pq(S) 
around E = 0.26 as reported in Fig. ^ Error bars are of 
the size of the circles. The System corresponds to A'^ = 230 
bosons and U — 280. See text for details. 



B. The band profile 

The fingerprints of classically chaotic dynamics can be 
found also in the band-structure of the perturbation ma- 
trix B. As we will show below the latter is related to the 
fluctuations of the classical motion. This is a major step 
towards a RMT modeling. 

Consider a given ergodic trajectory (I(t),ip(t)) on the 
energy surface 7Y(/(0), </3(0); fco) — E (with N = const.). 
We can associate with it a stochastic-like variable 



(16) 



which can be seen as a generalized force. For the BHH 
([5]) this is simply given by the perturbation term i.e. 



(17) 
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FIG. 6: (Color online) The power-spectrum of the classical 
trimer ^ at energy E = 0.26, U = 280, and Ao = 0.053. The 
classical cut-off frequency lo^i = iliaiU « 280 is indicated by 
vertical dashed lines. 



which corresponds to a momentum boost since it changes 
the kinetic energy [68j . It may have a non-zero average, 
i.e. a "conservative" part, but below we are interested 
only in its fluctuations. 

In order to characterize the fluctuations of !F{t) we 
introduce the autocorrelation function 



(18) 



where f = Ut is a rescaled time. The angular brackets 
denote an averaging which is either micro-canonical over 
some initial conditions (l(0),ip(0)) or temporal due to 
the assumed ergodicity. 

For generic chaotic systems (with smoothly varying po- 
tentials), the fluctuations are characterized by a short 
correlation time f^,!, after which the correlations are neg- 
ligible. In generic circumstances f„i is essentially the er- 
godic time. For our system we have found f^,, ~ 27r [see 

Eq. m]. 



The power spectrum of the fluctuations C(aj) is given 
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by a Fourier transform: 



C{f) e"^^df 



(19) 



and for the case of the trimer ^ is shown in Fig. [6] 
We see that C{ui) has a (continuous) frequency support 
which is bounded by w^i ~ 1 corresponding to iv^i fa 280 
(indicated by dashed vertical lines in Fig. |6|. The cut- 
off frequency lu^i is inverse proportional to the classical 
correlation time, i.e. 



27T 



(20) 



These characteristics of the power spectrum are universal 
for generic chaotic systems. Finally, we see that within 
the frequency support the power spectrum C{uj) is struc- 
tured, reflecting system-specific properties of the under- 
lying classical dynamics. 

The classical power spectrum C{uj) is associated with 
the quantum mechanical perturbation matrix B accord- 
ing to the following semiclassical relation [591 IZO] 



(|B„ 



— G uj 

U2tt 



En — E„ 



(21) 



Hence the matrix elements of the perturbation matrix B 
are extremely small outside a band of width 



b = huj,,/A « huj,,N/U 



(22) 



In the inset of Fig. [7] we show a snapshot of the per- 
turbation matrix |B„mp which clearly exhibits a band- 
structure. In the same figure we also display the scaled 
quantum band profile for N = 230. The agreement with 
the classical power spectrum C{uj) is excellent. We have 
checked that the relation is very robust [l9 l [29 1 154] 
and holds even for moderate number of bosons iV « 50. 
Combining Eqs. (12) and (22 1 with a)^,, w 1 (see above) 



and the definition of b we find for the chaotic regime 
around E ~ 0.26 that b ^ 0.6N which is confirmed by 
the numerics. 

It is important to realize that upon quantization we 
end up with two distinct energy scales ^T9\ l29l l54] . One is 
obviously the mean level spacing A ~ 1 /N (see Eq. ( 12 1 ) 
which is associated with the unperturbed Hamiltonian. 
The other energy scale is the bandwidth 



Ah 6A cx [/ 



(23) 



which contains information about the power spectrum 
of the chaotic motion and is encoded in the perturba- 
tion matrix B. The latter energy scale is also known in 
the corresponding literature as the "non-universal" en- 
ergy scale [TIj, or in the case of diffusive motion, as the 
Thouless energy [75]. One has to notice that deep in the 
semiclassical regime N ^ oo these two energy scales dif- 
fer enormously from one another. We shall see in the 
following sections that this scale separation has dramatic 
consequences on the theory of wavepacket dynamics. 
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FIG. 7: (Color online) The band profile {2nU/N^A) ■ |B„„p 
versus lo = {En — Em)/fi is compared with the classical power 
spectrum C{lj). The number of particles is N = 230 and 
Ao = 0.053. Inset: a snapshot of the perturbation matrix 



C. Distribution of matrix elements of the 
perturbation operator 

We further investigate the statistical properties of the 
matrix elements B„,„ of the perturbation matrix by 
studying their distribution. RMT assumes that upon 
appropriate "unfolding" they must be distributed in a 
Gaussian manner. The unfolding aims to remove system 
specific properties and to reveal the underlying universal- 
ity. It is carried out by normalizing the matrix elements 
with the local standard deviation a = \J (|B„,„p) related 
through Eq. (21 1 with the classical power spectrum C{lo). 



The existing literature is not conclusive about the dis- 
tribution of the normalized matrix elements w — B„m/tT. 
Specifically, Berry [73] and more recently Prosen [TOl |74] , 
claimed that V{w) should be Gaussian. On the other 
hand, Austin and Wilkinson [7S] have found that the 
Gaussian is approached only in the limit of high quan- 
tum numbers while for small numbers, i.e., low energies, 
a different distribution applies, namely 



r( 



V7rivr( 



2 I 



W 



(7V-3)/2 



(24) 



This is the distribution of the elements of an A^- 
dimensional vector, distributed randomly over the sur- 
face of an A-dimensional sphere of radius ^/N. For 
N oo this distribution approaches a Gaussian. 

In Fig. [s] we report the distribution V{w) for the el- 
ements of the perturbation matrix B. The dashed line 
corresponds to a Gaussian of unit variance while the cir- 
cles are obtained by fitting Eq. (|24| to the numerical data 
using iV as a fitting parameter. Although we are deep in 
the semiclassical regime (i.e. N — 230), none of the above 
predictions describes in a satisfactory way the numerical 
data. We attribute these deviations to the existence of 
small stability islands in the phase space. Trajectories 
started in those islands cannot reach the chaotic sea and 
vice versa. Quantum mechanically, the consequence of 
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FIG. 8: (Color online) Distribution of rescaled matrix ele- 
ments w around E = 0.26 rescaled with the averaged band 
profile. The dashed red line corresponds to the standard nor- 
mal distribution while the circles (o) correspond to a best fit 
from Eq. ( 24 1 with a fitting parameter A'' = 342. The system 
corresponds = 230, U = 280. 
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FIG. 9: (Color online) Schematic representation of the two 
notions of the kernel P{n\m). Left: Projection of one per- 
turbed eigenstate \n{ko + 5k)) (blue level) on the basis |m(fco)) 
of the unperturbed Hamiltonian. Averaging over several \n') 
states around energy En yields the averaged shape of eigen- 
functions (ASoE). Right: Alternatively, if P(n\m) is regarded 
as a projection of one unperturbed eigenstate \m) (blue level) 
on the basis |n) of the perturbed Hamiltonian and averaged 
over several states around Em , it leads to the local density 
of states (LDoS). 



this would be vanishing matrix elements B„m which rep- 
resent the classically forbidden transitions. 



D. RMT modeling 

More than 50 years ago, E. P. Wigner [62l |63] pro- 
posed a simplified model to study the statistical prop- 
erties of eigenvalues and eigenfunctions of complex sys- 
tems. It is known as the Wigner banded random ma- 
trix (WBRM) model. The corresponding Hamiltonian 
is given by Eq. ([s]) where B is a banded random matrix 
[751 [771175]. This approach is attractive both analytically 
and numerically. Analytical calculations are greatly sim- 
plified by the assumption that the off-diagonal terms can 
be treated as independent random numbers. Also from 
a numerical point of view it is quite a tough task to cal- 
culate the true matrix elements of B. It requires a pre- 
liminary step where Hq is diagonalized. Due to memory 
limitations one ends up with quite small matrices. For 
example, for the Bose-Hubbard Hamiltonian we were able 
to handle matrices of final size J\f = 30,000 maximum. 
This should be contrasted with RMT simulations, where 
using self-expanding algorithm [121 HZl [73] we were able 
to handle system sizes up to A/" = 1, 000, 000 along with 
significantly reduced CPU time. We would like to stress, 
however, that the underlying assumption of the WBRM, 
namely that the off-diagonal elements are uncorrelated 
random numbers, has to be treated with extreme care. 
The applicability of this model is therefore a matter of 
conjecture which we will test in the following sections. 

In fact, the WBRM model involves an additional sim- 
plification. Namely, one assumes that the perturbation 
matrix B has a rectangular band profile of bandwidth 



b. A simple inspection of the band profile of our BHH 
model (see Fig. [7| shows that this is not the case. We 
eliminate this simplification by introducing a RMT model 
that is even closer to the dynamical one. Specifically, 
we generate the matrix elements B„„^ from a Gaussian 
distribution with a variance that is given by the classi- 
cal power spectrum according to Eq. (21 1. Thus the 
band-structure is kept intact. This procedure leads to 
a random model that exhibits only universal properties 
but lacks any classical limit. We will refer to it as the 
improved random matrix theory model (IRMT). 



V. LOCAL DENSITY OF STATES AND 
QUANTUM-CLASSICAL CORRESPONDENCE 

As we change the parameter 5k in the Hamiltonian (|5| , 
the instantaneous eigenstates {|?T.(fc))} undergo struc- 
tural changes. Understanding these changes is a cru- 
cial step towards the analysis of wavepacket dynamics 
[2ni [5H . This leads to the introduction of the "kernel" 



P{n\m) = |(n(fco + (5fc)|m(fco))|^ 



(25) 



which can be interpreted in two ways as we schemat- 
ically depict in Fig. [9] If regarded as a function of 
TO, P(n\m) represents the overlap of a given perturbed 
eigenstate |n(A:o -I- 5k)) with the eigenstates |TO(fco)) of 
the unperturbed Hamiltonian. The averaged distribution 
P{r) is defined by r = n — m, and averaging over several 
states with roughly the same energy En yields the av- 
eraged shape of eigenfunctions (ASoE). Alternatively, if 
regarded as a function of n and averaging over several 
states around a given energy the kernel P{r) rep- 
resents up to some trivial scaling and shifting the local 
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FIG. 10: (Color online) The kernel P{n\m) of the BHH plot- 
ted as a function of the perturbed energies En (LDoS repre- 
sentation) and for various perturbation strengths Sk > (5feqm- 
The averaged shape of eigenfunctions is given by the same 
kernel P{n\m) and is obtained by just inverting the energy 
axis. Here, A*' = 70, and Ao = 0.053. 
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FIG. 11; (Color online) The quantal profile P{n\m) as a 
function of En — -Em' for the BHH model is compared with 
Pprt and with the corresponding Pirmt of the IRMT model. 
The perturbation strength Sk is in (a) standard perturbative 
regime Sk — 0.05 and (b) extended perturbative regime Sk = 
0.3. The system corresponds to N = 230, U = 280 and 



ko = 15. Here Skn 



: 0.09 and Sk^ 



1.02. 



density of states (LDoS): 

P{E\m) = Y,Mk)\m{ko))\^5{E - E.^) . (26) 



Its line-shape is fundamental for t he understanding of the 
associated dynamics (see Sec. VI I, since its Fourier trans- 

In 



form is the so-called "survival probability amplitude" 
the following we will focus on the LDoS scenario. 



A. Parametric Evolution of the LDoS 

An overview of the parametric evolution of the aver- 
aged P(n\m) is shown in Fig. 



delta function for 6k — 0, the profile P{n\m) starts to 
develop a non-perturbative core as Sk increases above 
some critical value Sk^„-^. For even stronger perturba- 
tions, P{n\m) spills over the entire bandwidth Af,. We 
will show that if Sk exceeds another critical value (5fcprt, 
the LDoS develops classical features. In the following we 
will identify the above parametric regimes and discuss 
the theory of P{n\m) in each one of them. 



1. The perturbative regimes 

We start with the discussion of the perturbative 
regimes. We distinguish between two cases: 

The Standard Perturbative Regime: The simplest case 
is obviously the first order perturbation theory (FOPT) 
regime where, for P(n\m), we can use the standard text- 



book approximation Pfopt('^|w) ^ 

Pfopt("-|to) 



i 1 for n = m, while 

Bo77,7j, I 



{En — Em)' 



(27) 



for n 7^ m. The border Sk^^„ for which Eq. (27 1 describes 
the LDoS kernel, can be found by the requirement that 
only nearest-neighbor levels are mixed by the perturba- 
tion. We get 



SK 



A/a cx 



U 
Ar3/2 



(28) 



10 Beginning as a 



where for the rhs. of Eq. ( 28 ) we have used the scaling re- 
lations for A and a (see Eqs. (12) and (21 1). In Fig. Hi 
we report our numerical results for the BHH, together 
with the perturbative profile Pj^oninlm) obtained from 
(27 1 and the outcome of the IRMT modeling. The 



FOPT Eq. (|27| has as an input the classical power spec- 

( |2l| can be used in order to 
All three curves fall on 



trum C(w) which via Eq. 
evaluate the band profile 
top of one another. 

Extended Perturbative Regime: If Sk > Sk^^ but not 
too large then we expect that several levels are mixed 
non-perturbatively. This leads to a distinction between 
a "core" of width F which contains most of the proba- 
bility and a tail region which is still described by FOPT. 
This non-trivial observation can be justified using per- 
turbation theory to infinite order. It turns out that the 
non-perturbative mixing on the small scale F of the core 
does not affect the long-range transitions [Ml [80] that 
dictate the tails. Therefore we can argue that a reason- 
able approximation is |54j 



Pprt(n|m) 



Sk^ \Br 



{En — En 



F2- 



(29) 



10 




FIG. 12: (Color online) Various measures of the spreading 
profile for the BHH and the IRMT model: the quantal spread- 
ing 5E (black line), the quantal spreading SEj^^t of the per- 
turbative profile given by Eq. ( |29[ l (red line), the spreading 
SEmMT obtained from the IRMT modeling, the analytical 
spreading 5iJana obtained from (34 1 (green line), the core- 
width r (orange line), and the classical spreading SE^i (blue 
line). The dashed line has slope one, while the dash-dotted 
line has slope two and are drawn to guide the eye. The sys- 
tems correspond to A'^ = 70 bosons, ko — 15 and U — 280. 
See text for details. 

Our numerical data, reported in Fig |ll| 3, indicate again 
an excellent agreement with the theoretical prediction 
( 29 1 . At the same time, we observe that also the proposed 
IRMT describes quite nicely the actual profile P{r) . Note 
that the resulting line-shape is strikingly difi^erent from a 
Wigner Lorentzian (as predicted by the traditional RMT 
modeling) and is rather governed by the semiclassical 
structures of the band profile |B„,„p. Instead, a Wigner 
Lorentzian would be obtained if the band profile of the 
perturbation matrix were flat. 

The core-width F is evaluated by imposing normaliza- 
tion on Pprt("-|'7i) |29| . Our numerically evaluated F is 
reported in Fig. [12 We see that for very small 6k we get 
that F ^ A. In this case, the expression (29 1 collapses to 
the FOPT expression (27 1. In fact, the inequality F < A 
can be used in order to estimate the limit 6k^^ of the 
validity of FOPT. As soon as we enter the extended per- 
turbative regime, we find (see Fig. 12 1 that F grows as 

Sk^YxA . (30) 



F oc 



The core-width F (and thus Eq. (|29| for the LDoS) is 
meaningful only as long as we have F < A^, i.e. as long 
as we can distinguish a core-tail structure. This condition 
allows us to evaluate the perturbative border dkp,^: 



u 

N 



(31) 



In our numerical analysis we have defined Sk^^ as the 
perturbation strength for which 50% of the probability 
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13: (Color online) The parameters (a) Skqm and (b) 
for various U, N, and for Aq = 0.053. A nice scaling in 



FIG. 

accordance with Eqs.(28l and (311 is observed 



remains at the original site but we have checked that the 
condition F = A gives the same result. For determining 
(5fcp,t we use the following numerical procedure: We cal- 
culate the spreading 5E = \/j2n -P(?i|"^)(£^m^ — En) 
P{r). Next we calculate 5E^^,{6k), using Eq.p9|)) 



2 of 
This 

quantity always saturates for large 5k because of having 
a finite bandwidth. We compare it to the exact 6E{6k) 
and define Jfcprt, for instance, as the 80% departure point. 



In Fig. 13 we present our numerical data for (5fcq„i and 



dkprt by making use of the scaling relations ( 28 1 and ( 31 1 . 
A nice overlap is evident, confirming the validity of the 
above expressions. 



2. The non-perturbative regime 

For Sk > Jfcpit the core spills over the bandwidth and 
therefore perturbation theory, even to infinite order, is 
inapplicable for evaluating P(n\m). In this regime, we 
observe that also the IRMT fails to reproduce the actual 
line-shape of P{n\m) as can be seen from Fig. 14 1. In 
fact, RMT modeling leads to a semicircle 



P{n\m) = l/{2TTA)^4-{{En~Em.)/Ay 



(32) 



as predicted by Wigner [55]. 

What is the physics behind the LDoS line-shape in 
the non-perturbative regime? Due to the strong per- 
turbations many levels are mixed and hence the quan- 
tum nature becomes "blurred" . Then, we can approxi- 
mate the spreading profile by the semiclassical expression 

[MllHollHT] 



Ps^{n\m) 



(33) 



where Pm(/, (p) and Pn{I,f) are the Wigner functions 
that correspond to the eigenstates |m(fco)) and \n{k)) 
respectively. In the strict classical limit p can be approx- 
imated by the corresponding micro-canonical distribu- 
tion p cx S(E—T-l({Ii}, {ipi})) determined by the energy 
surface E. The latter can be evaluated by projecting 
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FIG. 14: (Color online) Upper panel: The kernel P(n\m) 
(LDoS representation) in the non-perturbative regime Sk = 
10 for iV = 230 and A = 0.053. The results of the BHH model 
(solid black line) are compared with Pprt (red dashed line), 
PiBMT of the IRMT model (solid green line), and the classical 
profile Pel (blue line with o). Lower panel: A time series E{t) 
which leads to the classical profile Pci(-E) (see text for details). 



the dynamics generated hy Ho{{Ii},{fi}) = Eq onto the 
Haniiltonian H({/J, {(/jJ) = E{t). 

In Fig. [14]3 we plot the resulting E{t) = n{Iit), ip{t)) 
as a function of time for the Hamiltonian ^ . The clas- 



sical distribution Pe,(n|m) is constructed (Fig.[14ji) from 
E(t), by averaging over a sufficiently long time. The 
good agreement with the quantum profile P(n\m) is a 
manifestation of the detailed quantum-classical corre- 
spondence which affects the whole LDoS profile in the 
non-perturbative regime. 

Coming back to the failure of the IRMT approach, we 
are now able to understand it formally from the scaling 
relation ( [si] ) of the perturbative border Sk-p^^ ~ tj /N. 
Specifically, we observe that the non-perturbative limit 
can be approached either by increasing the perturbation 
strength Sk or, alternatively, by keeping Sk constant and 
increasing N. As we have seen before increasing A'^ means 
to approach the classical limit (keeping U — const.). On 
the other hand, it is clear that the IRMT model lacks a 
classical limit! Therefore, we cannot expect it to yield 
a correct description of P{n\m) in that regime. Instead, 
for 5k > (5fcprt the LDoS is completely dictated by semi- 
classical considerations as can be seen from Fig. [T4k. 



152] . The two types of QCC are defined as follows: (a) 
detailed QCC means P{r) ~ ^'ci('') while (b) restricted 
QCC means SE^^, « SE,^. 

Obviously restricted QCC is a trivial consequence of 
detailed QCC, but the converse is not true. It turns out 
that restricted QCC is much more robust than detailed 
QCC. In Fig. 12 we see that the dispersion SE^^ of ei- 
ther P{r) or -P[rmt(^) is almost indistinguishable from 
6E^i. In fact, this agreement of the second moment SE 



persists also for the case of the perturbative profile ( 29 1 . 
This is quite remarkable because the corresponding LDoS 
profiles (quantal, perturbative, IRMT and classical) can 
become very different! 

The possibility of having restricted QCC was pointed 
out in |54l I81j in the frame of quantum systems with 
chaotic classical limit. A simple proof presented in 
Ref. [54] indicated that the variance of P{r) is deter- 
mined by the first two moments of the Hamiltonian in 
the unperturbed basis i.e. 



6E' 



{m\H \m) 
Sk^ \(m\B^\m) 



lm\H\m) 



E 



{m\B\m) 



IB, 



(34) 



Having a 5E^^ that is determined only by the band pro- 
file, is the reason for restricted QCC, and is also the 
reason why restricted QCC is not sensitive to the RMT 
assumption. 



VI. WAVEPACKET DYNAMICS 

We now turn to the time-dependent scenario of the 
wavepacket dynamics which is related to the response of 
a system to a rectangular pulse. Its physical realization in 
the framework of the BHH has been described in Section 

En] 

In the next subsections we will discuss the time- 
evolving energy profile in each of the three (5fc-regimes 
which we have identified in the frame of the LDoS study. 
We sta rt our analysis with the classical dynamics (Sub- 
VI A[ ) and then turn to the evolution of the quan- 



section 

tum profile Pt(r) (Subsection VIB I. In the same subsec- 
tion we will present an analysis of the IRMT and semi- 
classical modeling and identify both their weakness and 
regimes of validity. 



A. Classical Dynamics 



B. Restricted vs. Detailed Quantum-Classical 
Correspondence 

It is important to distinguish between detailed and 
restricted quantum-classical correspondence (QCC) [THl 



The classical picture is quite clear: The initial prepara- 
tion is assumed to be a micro-canonical distribution that 
is supported by the energy surface Hq{I,(p) = E{0) — 
Eng where the Hamiltonian is given by Eq. Taking 
7i(A) to be a generator for the classical dynamics, the 
phase-space distribution spreads away from the initial 



12 




J I 1— 



FIG. 15: (Color online) The classical energy spreading 5Eai{t) 
for the BHH (normalized with respect to the perturbation 
strength Sk and the boson number N) is plotted as a function 
of time. The dashed line has slope one and is drawn to guide 
the eye. 



surface for t > 0. "Points" of the evolving distribution 
move upon the energy surfaces of 7i(/, </?). Thus, the en- 
ergy E{t) — Ho{I{t),ip{t)) of the evolving distribution 
spreads with time. We are interested in the distribution 
of E{t) of the evolving "points" . 

A quantitative description of the classical spreading is 
easily obtained from Hamilton's equations: 



dE (t) 
dt 



m 

Ik 



= -5kf{t)T{t) 



(35) 



5E,,{t) = Skx v/2[C(0) - C{t)] 



where [-Jpb indicates the Poisson Brackets and f{t) is a 
rectangular pulse i.e. f{t') = 1 for < t' < t. Inte- 
grating the previous expression and then taking a micro- 
canonical average over initial conditions we get for the en- 
ergy spreading the classical linear response theory (LRT) 
expression 

SEci:^^ ; t < Tci 
SEci ; t>Tci' 
(36) 

In the last step, we have expanded the correlation func- 
tion for t < Tci as C{t) « C(0) - ^C"{0)t^. For t > t^, 
due to ergodicity, a "steady-state distribution" appears, 
where the evolving "points" occupy an "energy shell" in 
phase-space. The thickness of this energy shell equals 
dE^,. Thus, the classical dynamics is fully characterized 
by the two classical parameters and SE^i. 

Figure [15] shows the scaled classical energy spreading 
6E^,{t)/{N Sk) for the BHH. The heavy dashed line has 
slope one and is drawn to guide the eye. In agreement 
with Eq. (36 1 we see that 6E^i{t) is first ballistic and then 
saturates at r^i « 2tt/U = 0.02. 

One can also calculate the entire classical evolving pro- 
file Pc\{t). Using a phase-space approach similarly to the 



onto Ho yields a set of energies {E}t=t whose distribu- 
tion [92] constitutes the spreading profile -Pci(^) at time 



t. We will discuss P^i{t) in Subsection VIC 



B. Quantum Dynamics 

Now we would like to explore the various dynamical 
scenarios that are generated by the Schrodinger equation 
for a„(t) — {n\tp{t)). Namely, we want to solve 



dan 
~dt 



-E„ 



(37) 



starting with an initial preparation a„ = S„m at t=0, i.e. 
an eigenstate of the unperturbed system. We describe 
the energy spreading profile for i > by the transition 
probabihty kernel Pt{n\m) — (|a„(t)p). The angular 
brackets stand for averaging over initial states (m) be- 
longing to the energy interval 0.25 < E^ < 0.27. We 
characterize the evolving distribution using the various 
measures introduced in subsection IIII Al If the evolution 
is classical-like then -according to the classical analysis 
presented previously- Pt(n|TO) will be characterized by a 
single energy scale SE(t), meaning that any other mea- 
sure like SEcorcit) reduces (up to a numerical factor) to 
5E(t). We will use this criterion in the following in order 
to identify for which (Sfc-regimes the evolution is classical- 
like and for which ones it develops quantum features. 

An overview of the spreading profiles for three repre- 
sentative (5fc-strengths is given in Fig. [16] A qualitative 
difference in the spreading is evident: In Fig. [T6[ i the 
probability is mainly concentrated in the initial level for 
all times (standard perturbative regime). In Fig. 



16 



3 one 



can distinguish two different components in the Pt{n\ra), 
the "core" (characterized by 5Ecorc{t)) and the "tail" 
component (characterized by SE{t)), both of them be- 
ing smaller than the bandwidth (extended perturbative 
regime). For even stronger perturbations, the core spills 
all over the bandwidth (see Fig. [l6] :) and the dynamics 
is non-perturbative. In the following we discuss each of 
these regimes separately. 



1. The perturbative regimes 

For small perturbations 6k < dk^^ (see Fig. [16^ ) the 
probability is mainly concentrated in the initial level dur- 
ing the entire evolution. This is the FOPT (standard 
perturbative) regime where the perturbation mixes only 
nearby levels and little probability escapes to the tails. 

As the perturbation strength is increased (5fcq,„ < Sk < 



(Fig. 16 



LDoS case in Subsection V A2 we propagate up to time 
t under the Hamiltonian H, a large set of trajectories 
{E}t^Q that originally are supported by the energy sur- 
face 'Hq{I, (p) = E{t = 0) = Eno- Projecting them back in Fig. 17 1 together with the classical spreading SE^i{t) 



3 j) , levels within the bandwidth are mixed 
and one can distinguish two different components in the 
profile Pt(r): The core characterized by 6E^^^^{t), where 
most of the probability is concentrated, and the tail com- 
ponent, characterized by 5E(t). The latter is reported 
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b) 
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FIG. 16: (Color online) The profile Pt(r) of the BHH plotted 
as a function of time for various perturbation strengths Sk < 
5fcq„ (a), 5fcqm < Sk < (Jfcprt (b), 5k > Skp^t (c). Note the 
diflerent scale in (c). Here, N = 70, E = 0.26 and Aq = 0.053. 



The remarkable fact is that, as far as 5E{t) is concerned, 
the agreement with the classical result is perfect. This 
might lead to the wrong impression that the classical and 
quantum spreading are of the same nature. However, this 
is definitely not the case. 

In order to reveal the different nature of the quan- 
tum spreading in the perturbative regime we turn to the 
analysis of the core-width 5E^^^^{t) (see Fig. 17d). If the 



spreading were of classical type this would imply that 
the evolving profile would be characterized by a single 
energy scale, and thus 5E{t) ^ 5E^^^^{t). However, as 
can be seen in Fig. [T7)d this is certainly not the case: 
For 5k < Sk^,„ we have that SEcoreit) = A for all times 
while for 5fc,„ < Sk < dkp.t the core- width fulfills the 
inequalities A < SEcoi-dt) < 5E{t) < Af,. In fact, this 
separation of energy scales allows us to use perturbation 
theory in order to evaluate theoretically the evolving sec- 
ond moment of the energy distribution. We get for the 
transition probability from an initial state m to any other 
state n m 



Sk^ 



Pt{n\m) - — 2-|-"nm| 



(38) 



Here Ft{uj) = (ti;t)^-sinc^(a;i/2) is the spectral content 
of a constant perturbation of duration t, and sinc(a;) = 
sin(j:)/2:. Substituting the above expression in Eq. ([9| 
we get the LRT expression (36) for SE{t). We have 



also calculated the second moment resulting from the 
IRMT modeling. The outcome is reported in the in- 
set of Fig. [TT^ and shows that within the perturbative 
regime the IRMT modeling provides the same results (as 
far as the second moment is concerned) as the LRT cal- 
culations. Therefore we conclude that for Sk < Skp,-t 
the IRMT modeling, the LRT results, the classical re- 
sults SEci, and the quantum calculations for the second 
moment SE{t) of the BHH match one another. 

Encouraged by this success of LRT and the IRMT 
modeling to describe the second moment SEit) of the en- 
ergy spreading, we can further use them to evaluate the 
survival probability 'P{t). Assuming a Markovian picture 
of the dynamics, LRT predicts P2] 



V{t) = exp 



-Sk X / —C{uj)- 



2tt 



(39) 



which after substituting the spectral-content Ft{uj), can 
be re-written in the following form 



Vit) = exp 



°° duj ~, . , ^ f Lot 
— C{u) t tsmc^ - 



(40) 

For short times (t <^ t^,) during which the spreading is 
ballistic-like, the term tsmc^{iL;t/2) is broad compared to 
the band profile and can be approximated by t leading 
to 



r{t) = exp 



-C(r=0) X 



Skt 
IT 



(41) 
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FIG. 17: (Color online) ie/t panel: The (normalized) energy spreading SE{t) for the BHH and the IRMT model (inset) for three 
different perturbation strengths 5fc = 0.05<(5fcq„, (solid black line), 5k^^<Sk — 0.3 < Skp^t (dashed red line), and 5k — 5 > Sk^^t 
(dash-dotted green line). The classical expectation SEci{t) is represented in all three plots by a dashed blue line for comparison. 
In the inset the black dash-dotted lines have slope one and one-half respectively and are drawn to guide the eye. While for the 
BHH model one observes restricted quantum-classical correspondence in all regimes this is not the case for the IRMT model 
(inset): For perturbations Sk > Sk^^t the energy spreading 5E(t) exhibits a premature crossover to diffusive behavior. Right 
panel: The evolution of the corresponding core width 5Ecarc{t) for the BHH model is plotted. In the perturbative regimes 
one observes a separation of scales SEcai^it) < SE{t) < Ab, which is lost for strong perturbations Sk > Skp^t, where 5Ecaic{t) 
approaches more and more the classical expectation SEci{t). Here, A'' = 230, U — 280, E = 0.26 and Ao ~ 0.053. 



For longer times (t 3> r^i) on the other hand, the term 
tsinc^(ci;t/2) is extremely narrow and can be approxi- 
mated by a delta function S{uj). This results in a Fermi- 
Golden-Rule (FGR) decay 



V{t) = exp 



C{uj = 0) X i 



(42) 



which can be trusted as long a,s 'P{t) ~ 1. This can be 

converted into an inequality t < iprt = ^TAf^) '^ci- 

In Fig. [18^ we plot our numerical results for the 
trimeric BHH model together with the theoretical ex- 



pectation ( 39 1 (we note that the outcome of the IRMT 



modeling matches exactly the results of the LRT and 
thus we do not overplot them). In both perturbative 
regimes we observe a short initial Gaussian decay (as im- 
plied by Eq. (41 1) which is followed by the exponential 
FGR decay. In the FOPT regime (inset of Fig. the 
entire decay until saturation is described by LRT. In the 



extended perturbative regime (see Fig. 18 1) the overall 



agreement is still pretty good. However, here the pertur- 
bative break time ip^t is shorter and one finds a deviation 



around the time t„ 



0.01. 



2. The non-perturbative regime 

Once we enter the non-perturbative regime Sk > ^fcpjt 
(see Fig. 16:), the core spills over the bandwidth and the 
separation of energy scales is lost, leading to 6E{t) ~ 
(5i?oorc(i) > Af, (see Fig. [TtId for Sk = 5). In this case 



the evolving energy distribution becomes totally non- 
perturbative. Still, for short times iprt = (^Tif^) "^ci < TcI, 
defined by the requirement that 'P{t) ~ 1 (see Eq. 
(41 1), the evolving probability kernel Pt{n\m) (and there- 



fore the spreading SE(t)) is described accurately by the 
FOPT expression ([38]). 

The remarkable fact is that although for t > tprt the 
evolving profile P(n|rn) is totally non-perturbative, this 
crossover is not reflected in the variance (see Fig. 17 1). 
The agreement with the LRT results of Eq.(|36| is still 
perfect. Instead, the crossover can be detected by study- 
ing other moments like SEcoroit) which acquire classi- 
cal characteristics, i.e. SEcore{t) ~ 6E{t) = SEc\{t) (see 
Fig. 17 d). Thus we are led to the conclusion [T5J [55] 



that the LRT predictions are not applicable while de- 
tailed QCC would possibly validate semiclassical consid- 
erations. We will examine this assumption more carefully 
in Subsection IVI CI 

What about the IRMT modeling? In the inset of 
Fig. |17^ we report the numerical results for the spread- 
ing 5E{t) of the IRMT model. We observe that as soon 
as we enter the non-perturbative regime, the spreading 
SE(t) shows a qualitatively different behavior than the 
dynamical BHH model. Namely, after an initial bal- 
listic spreading (taking place for times t < tprt), we 
observe a premature crossover to a diffusive behavior 
5E{t) — \/2DEt. The following heuristic picture can 
explain the diffusive behavior of the IRMT modeling. 
At < ~ tprt ^ Tci, the evolving distribution becomes as 
wide as the bandwidth, and we have 5E^^^^ ~ 5E ~ A;, 
rather than SE^„^^ <C SE <C A{,. Once the mechanism for 
ballistic- like spreading disappears, a stochastic-like be- 
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FIG. 18: (Color online) The survival probability for the BHH and three different perturbation strengths a)5k = 0.05 < (5A;q„ 
(inset), 5fcqm < 5k = 0.3 < Sk^rt (main figure), and h) 5k = 5 > 5kprt ■ The solid black line represent the exact numerical result 
while the dash-dotted red line is the LRT result ( 39 1 calculated using the IRMT model. The inset of panel (a) represents the 

0.1 (see 



FOPT regime, while the main figure corresponds to the extended perturbative regime. Here the break time is t^^t 
Subsection (VIB2I. In the non-perturbative regime (b), the LRT breaks down close to the calculated break time 
In this panel we superimpose the Fourier transform of the LDoS as blue circles. 
N = 230, U = 280, E = 0.26 and Ao = 0.053. 



0.001. 

The agreement with 'P(t) is excellent. Here, 



havior takes its place. This is similar to a random- walk 
process where the step size is of the order A^, with tran- 
sient time tprt. 

The same deviations are observed for other observables 
as well. In Fig. \18\ ) we report our results for the survival 
probability in the non-perturbative regime. We find that 
the IRMT modeling (which for short times gives the same 
results as LRT-not shown in the figure as they are indis- 
tinguishable from the IRMT results) breaks down after 
an initial Gaussian decay ( 41 1 which holds up to a break 
time iprt ~ 0.001. Instead, the behavior of 7'(i) can be 
obtained by a Fourier transform of the LDoS. Specifically, 
we have that 



FOPT -or equivalently 



V{t) ^ \{n{ko)\e-'"^'^'^''\n{ko)) 



|(m(fc)|n(fco))P 

2 



(43) 



where P{E\m) is given by Eq. ( 26 1 . In Fig. 18 d we super- 
impose the outcome of Eq. ( 43 ((see blue circles) together 
with the survival probability evaluated by the numerical 
integration of the Schrodinger equation. An excellent 
agreement is evident. 



C. Detailed versus restricted QCC 

In the previous subsection we have assumed that 
the evolving wavepacket is developing detailed QCC in 
the non-perturbative regime and for times t > tpi-t = 



-jjP j (for earlier times 

IRMT considerations- apply). 

In Fig. 19 we report four snapshots of the evolving 
quantum mechanical profile (black lines). In the same 
figure we report the IRMT results (red lines) together 
with the classical calculations (blue lines with o). As we 
have discussed above we distinguish two phases in the 
evolution: For t < Tc\ the IRMT modeling (or equiv- 
alently the FOPT) is applicable while for t > Tc\ the 
evolving profile is described by its classical counterpart 
Pci(t). During this second phase, the evolution predicted 
by the IRMT is diffusive leading to a Gaussian shape for 
Pt{n\m). 



VII. CONCLUSIONS 

In this paper we have studied the evolving energy dis- 
tribution of a three-site ring-shaped Bose-Hubbard model 
in the chaotic regime. The evolution is triggered by a 
change dk in the tunneling rate k between neighboring 
lattice sites which in the context of ultra-cold atoms in 
optical lattices is realized by a change in the intensity 
of the trapping laser field. The specific scenario that 
we have analyzed in detail is the so-called wavepacket 
dynamics in energy space corresponding to a constant 
driving pulse of finite duration t. 

We followed a three-fold approach to the problem 
which combines purely quantum mechanical as well as 
semiclassical and random matrix theory considerations. 
This enabled us to identify both the strengths and limi- 
tations of each method. 

We find the appearance of three dynamical 5fe-regimes: 
The standard perturbative {5k < 5fcq„, oc U /N'^/'^), the 
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FIG. 19: (Color online) Snapshots of the evolving quantum profile Pt{r) obtained from the BHH (black line) and the IRMT 
model (red line) as well as the classical profile Pt^r) (blue line with o) in the non-perturbative regime 5k = 5 > Skp^t plotted 
against the energy difference E — Eq. After the quantal transition period t ~ 0.002 (see Fig. |17| 3) there is no scale separation 
between the core and the tail component and one observes overall detailed QCC. However, the initially excited component |no) 
decays slower in the quantum case. Here, A*' = 230, U = 280, E = 0.26 and Ao = 0.053. 



extended perturbative ((5fcq„, < 6k < (Jfcpjt oc U/N) and 
the non-perturbative regime {6 > Skp,^). The first two 
regimes can be addressed using LRT or RMT calcula- 
tions. In contrast, the last regime requires a combination 
of LRT/RMT calculations and semiclassical considera- 
tions. The former approach describes the evolving energy 
distribution for short times while the latter applies for 
longer times. Interestingly enough we have found that 
the variance 5E'^{t) of the evolving energy distribution 
shows a robust quantum-classical correspondence for all 
(5A;-values, while other moments exhibit this QCC only in 
the non-perturbative regime identified with the classical 
limit. In this regime, even an improved RMT modeling 
fails to describe the long time behavior of 6E{t) lead- 
ing to a premature crossover from ballistic to diffusive 
behavior. 

The motivation of the present study is driven both by 
theoretical and experimental considerations. On the fun- 
damental level, we would like to understand the manifes- 
tation of quantum-classical correspondence in the con- 



text of quantum chaotic dynamics, where chaos enters 
not due to geometrical considerations ("chaotic" shape 
of the trap) but due to many-body interactions [83^. At 
the same time, our results are also of immediate rele- 
vance to various branches of physics. For example, in the 
framework of ultra-cold atoms loaded in optical traps one 
is interested in understanding measurements of the en- 
ergy absorption rates induced by potential modulations 
[23111; 22, 23 ,121 US US]- Another application arises in 
molecular physics: As mentioned in Section |ll] the Bose- 
Hubbard Hamiltonian also models bond-excitations in 
small molecules [571 [53]. In this respect, the wavepacket 
dynamics investigated here describes the vibrational en- 
ergy redistribution of an initial excitation [84 . 

As far as the experimental realization of our study is 
concerned, micro-traps [5] are promising candidates for 
such time-dependent potentials [85] while optical lattices 
have already been successfully used in similar setups. 
Specifically, the studied dynamical scenario is readily im- 
plemented by changing the intensity of the laser field us- 
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ing a simplified version of the experiments of the Zurich 
group ^20l [21]. In contrast to the periodic modulation 
presented there, the optical lattice depth has to be altered 
in a step-like manner. Such experiments have been suc- 
cessfully performed by Greiner et al. [86j where the inten- 
sity of the trapping laser field was suddenly raised. The 
raise time was achieved to be much faster than the tun- 
neling time between neighboring sites but slow enough 
as not to excite higher vibrational modes of the wells. 

Concerning the measurement of the energy distribu- 
tion Pt (E) and the associated absorption of energy due 
to the driving, various techniques may be applied. Using 
standard time-of-flight measurements one can determine, 
for example, the release energy of the condensate and the 
momentum distribution of the atomic cloud which we ex- 
pect to provide the relevant information on the variance 
5E'^{t) of the energy distribution. Another possibility 
is to probe the Pt{E) via phase diffusion measurements 
[86j . Experimentally, the BEC can be prepared (almost) 
in one eigenstate. The driving pulse induces a broadening 



in the energy distribution leading to (decaying) oscilla- 
tions in the contrast {b\bi+i+bih\j^i) between neighboring 
sites. We expect that the functional form of the decay 
can be directly related to the core width V and thus be 
used to detect the three parametric (5fc-regimes. While 
these measurements are in principle sensitive to decoher- 
ence due to residual interaction with the non-condensed 
atoms, we note here that for two-site systems coherence 
times of several hundred milliseconds were observed [93] 
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